Fourier re-binning of time-of-flight positron emission tomography data

ABSTRACT

Fast reconstruction methods are provided for 3D time-of-flight (TOF) positron emission tomography (PET), based on 2D data re-binning. Starting from pre-corrected 3D TOF data, a re-binning algorithm estimates for each transaxial slice the 2D TOF sinogram. The re-binned sinograms can then be reconstructed using any algorithm for 2D TOF reconstruction. A TOF-FORE (Fourier re-binning of TOF data) algorithm is provided as an approximate re-binning algorithm obtained by extending the Fourier re-binning method for non-TOF data. In addition, two partial differential equations are identified that must be satisfied by consistent 3D TOF data, and are used to derive exact re-binning algorithms and to characterize the degree of the approximation in TOF-FORE. Numerical simulations demonstrate that TOF-FORE is more accurate than two different TOF extensions of the single-slice re-binning method, and suggest that TOF-FORE will be a valuable tool for practical TOF PET in the range of axial apertures and time resolutions typical of current scanners.

CLAIM OF PRIORITY FROM RELATED APPLICATION

This application claims priority under 35 U.S.C. § 119(e) from copendingProvisional Application Ser. No. 60/638,532, filed Dec. 22, 2004.

FIELD OF THE INVENTION

The present invention generally relates to nuclear medicine, and systemsfor obtaining nuclear medicine images. In particular, the presentinvention relates to systems and methods for reconstructing nuclearmedicine images from time-of-flight (TOF) positron emission tomography(PET) data.

BACKGROUND OF THE INVENTION

Nuclear medicine is a unique medical specialty wherein radiation is usedto acquire images which show the function and anatomy of organs, bonesor tissues of the body. Radiopharmaceuticals are introduced into thebody, either by injection or ingestion, and are attracted to specificorgans, bones or tissues of interest. Such radiopharmaceuticals producegamma photon emissions which emanate from the body and are captured by ascintillation crystal, with which the photons interact to produceflashes of light or “events.” Events are detected by an array ofphotodetectors, such as photomultiplier tubes, and their spatiallocations or positions are calculated and stored. In this way, an imageof the organ or tissue under study is created from detection of thedistribution of the radioisotopes in the body.

One particular nuclear medicine imaging technique is known as PositronEmission Tomography, or PET. PET is used to produce images fordiagnosing the biochemistry or physiology of a specific organ, tumor orother metabolically active site. Measurement of the tissue concentrationof a positron emitting radionuclide is based on coincidence detection ofthe two gamma photons arising from positron annihilation. When apositron is annihilated by an electron, two 511 keV gamma photons aresimultaneously produced and travel in approximately opposite directions.Gamma photons produced by an annihilation event can be detected by apair of oppositely disposed radiation detectors capable of producing asignal in response to the interaction of the gamma photons with ascintillation crystal. Annihilation events are typically identified by atime coincidence between the detection of the two 511 keV gamma photonsin the two oppositely disposed detectors, i.e., the gamma photonemissions are detected virtually simultaneously by each detector. Whentwo oppositely disposed gamma photons each strike an oppositely disposeddetector to produce a time coincidence event, they also identify a lineof response, or LOR, along which the annihilation event has occurred. Anexample of a PET method and apparatus is described in U.S. Pat. No.6,858,847, which patent is incorporated herein by reference in itsentirety.

After being sorted into parallel projections, the LORs defined by thecoincidence events are used to reconstruct a three-dimensionaldistribution of the positron-emitting radionuclide within the patient.In two-dimensional PET, each 2D transverse section or “slice” of theradionuclide distribution is reconstructed independently of adjacentsections. In fully three-dimensional PET, the data are sorted into setsof LORs, where each set is parallel to a particular detector angle, andtherefore represents a two dimensional parallel projection p(s, φ) ofthe three dimensional radionuclide distribution within the patient,where s corresponds to the distance of the imaging plane perpendicularto the scanner axis and φ corresponds to the angle of the detector planewith respect to the x axis in (x, y) coordinate space (in other words, φcorresponds to a particular LOR direction). Coincidence events areintegrated or collected for each LOR and stored as a sinogram. In thisformat, a single fixed point in f (x,y) traces a sinusoid in thesinogram. In each sinogram, there is one row containing the LORs for aparticular azimuthal angle φ; each such row corresponds to aone-dimensional parallel projection of the tracer distribution at adifferent coordinate along the scanner axis. This is shown conceptuallyin FIG. 1.

An event is registered if both crystals detect an annihilation photonwithin a coincidence time window τ (e.g., on the order of 4-5 ns),depending on the timing properties of the scintillator and the field ofview. A pair of detectors is sensitive only to coincidence eventsoccurring in the volume between the two detectors, thereby eliminatingthe need for physical collimation, and thus significantly increasingsensitivity. Accurate corrections can be made for the self-absorption ofphotons within the patient (i.e., attenuation correction) so thataccurate measurements of tracer concentration can be made.

The number of time coincidences detected per second within a field ofview (FOV) of a detector is the count rate of the detector. The countrate at each of two oppositely disposed detectors, A and B, can bereferred to as singles counts, or singles, S_(A) and S_(B). The timerequired for a gamma photon to travel from its point of origin to apoint of detection is referred to as the time of flight, or TOF, of thegamma photon. TOF is dependent upon the speed of light c and thedistance traveled. A time coincidence, or coincidence event, isidentified if the time difference between the arrival of signals in apair of oppositely disposed detectors is less than a coincidence timewindow τ.

As illustrated in FIG. 1, if an annihilation event occurs at themidpoint of a LOR, the TOF of the gamma photon detected in detector A(T_(A)) is equal to the TOF of the gamma photon detected in detector B(T_(B)). If an annihilation event occurs at a distance Δx from themidpoint of the LOR, the difference between T_(A) and T_(B) is Δt=2Δx/c,where c is the speed of light. If d is the distance between thedetectors, the TOF difference Δt could take any value from −d/c to +d/c,depending on the location of the annihilation event.

Time-of-flight (TOF) positron emission tomography (PET) (“TOF-PET”) isbased on the measurement of the difference Δt between the detectiontimes of the two gamma photons arising from the positron annihilationevent. This measurement allows the annihilation event to be localizedalong the LOR with a resolution of about 75-120 mm FWHM, assuming a timeresolution of 500-800 ps (picoseconds). Though less accurate than thespatial resolution of the scanner, this approximate localization iseffective in reducing the random coincidence rate and in improving boththe stability of the reconstruction and the signal-to-noise ratio (SNR),especially when imaging large objects.

TOF scanners developed in the early 1980s were used for research andclinical applications, but the SNR gain provided by the TOF measurementsof about 500 ps resolution was offset by the low stopping power of theBaF₂ and CsF scintillation crystals used in such scanners. Imagereconstruction for complete 2D TOF-PET data has been disclosed in theprior art. Early TOF scanners used a back project-then-filter (BPF)algorithm. The maximum-likelihood estimation algorithm (MLEM) also wasadapted for list-mode TOF data, and shown to provide improved imagequality compared to BPF. The increased computation required to processthe TOF data also led to the proposal of a faster iterative algorithmthat was directly applied to the back-projected TOF data.

In contrast with 2D TOF reconstruction, few studies have been devoted tothe 3D case, probably because the rise of 3D PET in the late 1980sovershadowed the interest for TOF. The Colsher filter and the 3Dre-projection algorithm for axially truncated data generalized the 3Dback projection filtered methods to TOF data.

However, the spatial resolution and sensitivity offered by those TOFsystems could not compete with the values achieved with BGO scanners. Asa result, TOF-PET almost completely disappeared from the scene in the1990s. Today, faster electronics and crystals such as LSO and LaBr₃reopen the prospect of exploiting the TOF information withoutcompromising other parameters such as the count rate, the sensitivity,and the energy and spatial resolutions. This prospect motivates thepresent invention of fast reconstruction strategies for 3D TOF-PET.

However, while extending MLEM or MAP (Maximum A posteriori Probabilitry)algorithms to 3D TOF-PET is conceptually straightforward, thecomputational load is an issue, because the number of data bins is nowequal to the number of LORs (which exceeds 10⁸ in 3D PET even afteraxial undersampling using the ‘span’ concept) multiplied by the numberof sampled TOF bins. What is needed is a feasible approach forreconstruction of 3D TOF-PET data.

SUMMARY OF THE INVENTION

The present invention solves the existing need in the art by providing ahybrid reconstruction method that re-bins the data onto a lowerdimensional space and then applies a conventional 2D reconstruction.This approach requires access to either pre-corrected data or correctionfactors (e.g., sensitivity, random coincidences, attenuation andscatter). The hybrid method offers a good compromise between imagequality and computational efficiency. The invention further makes thehypothesis that hybrid methods will also be effective for whole-bodyclinical TOF-PET imaging.

The first component of the hybrid approach according to the invention isthe re-binning algorithm. Starting from pre-corrected TOF projectiondata acquired by a volume (i.e., 3D) PET scanner, a re-binning algorithmestimates for each transaxial slice the ordinary 2D data set. There-binning algorithms according to the invention apply this procedureseparately for each TOF bin, and generate for each transaxial slice aTOF data set, which can be reconstructed using any analytic or iterativealgorithm for 2D TOF data.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will now be more fully described by way of example withreference to the accompanying drawings in which:

FIG. 1 is a diagram illustrating the relationship between PET projectiondata and a sinogram;

FIG. 2 is a diagram illustrating the concept of time of flight in PETimaging;

FIG. 3 is a composite diagram illustrating the parameterization ofTOF-PET data in accordance with the invention;

FIG. 4 is a graph illustrating normalized RMS error for re-binnedsinograms of an ellipsoid phantom as a function of slice index, for anumber of re-binning algorithms;

FIG. 5 is a diagram illustrating a cross-section of a cylindricalphantom as used in simulations in accordance with the invention;

FIG. 6 is a graph illustrating normalized RMS error for re-binnedsinograms of the cylindrical phantom at a FWHM of 120 mm;

FIG. 7 is a graph illustrating normalized RMS error for re-binnedsinograms of the cylindrical phantom at a FWHM of 75 mm;

FIGS. 8( a)-8(d) show reconstructed re-binned sinograms of a cylindricalphantom, according to various re-binning algorithms; and

FIGS. 9( a)-9(d) show axial sections of a stack of re-binned sinogramsof a cylindrical phantom, according to various re-binning algorithms.

DETAILED DESCRIPTION OF THE INVENTION

The present invention will now be described and disclosed in greaterdetail. It is to be understood, however, that the disclosed embodimentsare merely exemplary of the invention and that the invention may beembodied in various and alternative forms. Therefore, specificstructural and functional details disclosed herein are not to beinterpreted as limiting the scope of the claims, but are merely providedas an example to teach one having ordinary skill in the art to make anduse the invention.

Initially, the measured TOF data is parameterized as

$\begin{matrix}{{p_{t}^{m}\left( {s,\phi,z,\delta} \right)} = {\sqrt{1 + \delta^{2}}{\int_{- \infty}^{\infty}\mspace{7mu}{{\mathbb{d}{{lf}\left( {{{s\;\cos\;\phi} - {l\;\sin\;\phi}},{{s\;\sin\;\phi} + {l\;\cos\;\phi}},{z + {l\;\delta}}} \right)}}{h\left( {t,{l\sqrt{1 + \delta^{2}}}} \right)}}}}} & (1)\end{matrix}$where s and φ are the usual transaxial sinogram coordinates, z is theaxial coordinate of the mid-point of the LOR and δ=tan θ is the tangentof the angle θ between the LOR and a transaxial plane. These parametersare illustrated in FIG. 3. The range of these variables is the same asfor a non-TOF scanner.

For instance, for a cylindrical scanner with radius R_(d) and axialfield of view zε[0,L], the range is |s|≦R_(FOV), φε[0, π), and

$\begin{matrix}{{{\delta }\sqrt{R_{d}^{2} - s^{2}}} \leqslant z \leqslant {L - {{\delta }\sqrt{R_{d}^{2} - s^{2}}} - \frac{L}{2\sqrt{R_{d}^{2} - s^{2}}}} \leqslant \delta \leqslant \frac{L}{2\sqrt{R_{d}^{2} - s^{2}}}} & (2)\end{matrix}$where R_(FOV) denotes the radius of the cylindrical support of f. Notethat the integration variable l in equation (1) is the path lengthprojected onto the transaxial plane, and is related to the path length ralong the oblique LOR by r=l√{square root over (1+δ²)}=l/cos θ. For afixed pair (z, δ), the function p_(t) ^(m)(s, φ, z, δ), seen as afunction of s and φ, is referred to as an oblique sinogram.

In equation (1) the subscript t denotes the TOF bin, corresponding to asensitivity profile h(t, r) centered at position r=t along the LOR. TheTOF parameter is related to the difference Δτ between the arrival timesof the two photons by t=cΔτ/2, where c is the speed of light. There-binning algorithms according to the invention are derived formodified data defined by

$\begin{matrix}{{p_{t}\left( {s,\phi,z,\delta} \right)} = {\int_{- \infty}^{\infty}\mspace{7mu}{{\mathbb{d}{{lf}\left( {{{s\;\cos\;\phi} - {l\;\sin\;\phi}},{{s\;\sin\;\phi} + {l\;\cos\;\phi}},{z + {l\;\delta}}} \right)}}{{h\left( {t,l} \right)}.}}}} & (3)\end{matrix}$This parameterization is such that oblique sinograms with the same tvalue have the same TOF profile projected onto the transaxial plane, aproperty which facilitates re-binning. The simplest way to obtain themodified data is to make the approximationh(t,l)≅√{square root over (1+δ²)}h(t,l√{square root over (1+δ²)}),leading top _(t)(s,φ,z,δ)≅p _(t) ^(m)(s,φ,z,δ)  (4)or, alternatively, the better approximationh(t,l)≅√{square root over (1+δ²)}h(t√{square root over (1+δ²)},l√{squareroot over (1+δ²)}),leading top _(t)(s,φ,z,δ)≅p _(t) ^(m)√{square root over (1+δ2)}(s, φ, z, δ)  (5)As will be seen from the numerical results presented below, theseapproximations are accurate when the axial aperture δ is small (e.g.,smaller than 15°) and the TOF profile is sufficiently wide (e.g., a FWHMlarger than 500 ps), as for typical applications with current scanners.When the TOF profile is too narrow, or the axial aperture is too large,an exact expression can be used to calculate the modified data.

The SSRB (Single Slice Re-Binning) and TOF-SSRB equations below handleeach TOF bin separately and can be applied with an arbitrary TOF profileh(t, l). The Fourier re-binning algorithm presented below requiresmerging two opposite TOF bins t and −t. This leads to a symmetryrequirement h(t, l)=h(−t, −l), which is normally satisfied in practice.Finally, the exact re-binning equations disclosed below are valid onlyfor a shift invariant Gaussian profile.

The aim of a re-binning algorithm is to estimate, for each time bin t,the 2D TOF sinogram of each transaxial slice z within the axialfield-of-view of the scanner. This 2D TOF sinogram is defined byp _(reb,t)(s,φ,z)=p _(t)(s,φ,z, 0).  (6)

The re-binned sinograms therefore can be obtained directly from equation(6), by extracting the subset of the 3D data that corresponds to δ=0(or, in practice, to |δ| smaller than some threshold). A usefulre-binning algorithm, however, must incorporate the oblique sinogramsfor all available values of δ so as to optimize the SNR. The simplestre-binning algorithm is the straightforward extension of thesingle-slice re-binning (SSRB) algorithm for non-TOF data, and is basedon the approximationp _(reb,t)(s,φ,z)≅p _(t)(s,φ,z,δ).  (7)This relation is derived by neglecting the term proportional to δ in thethird argument of f in equation (3). An alternative is to take advantageof the approximate TOF localization, which allows a more accurateassignment of the LOR to a specific transaxial slice. Thus, using thefact that the TOF profile is maximum at l=t, the variable l in the thirdargument of f in equation (3) can to a first approximation be replacedby t, yielding the TOF-SSRB equation,p _(reb,t)(s,φ,z+tδ)=p _(t)(s,φ,z,δ).  (8)

Contrary to equation (7), the TOF-SSRB equation (8) becomes exact whenthe width of the TOF profile tends to zero (in which case, however,re-binning and reconstruction are no longer needed). Once a re-binningequation has been selected, a re-binned sinogram is obtained byaveraging the independent estimates provided by the redundant 3D data.For instance, if we use equation (7),

$\begin{matrix}{{p_{{reb},t}\left( {s,\phi,z} \right)} \simeq {\frac{1}{2{\delta_{\max}(z)}}{\int_{- {\delta_{\max}{(z)}}}^{\delta_{\max}{(z)}}\;{{\mathbb{d}\delta}\;{p_{t}\left( {s,\phi,z,\delta} \right)}}}}} & (9)\end{matrix}$where δ_(max)(z) determines the range of available oblique sinograms forslice z (see equation (2)). The expression is similar for the TOF-SSRBalgorithm.

The Fourier re-binning algorithm for non-TOF data was originally derivedby applying the stationary phase approximation to the 2D Fouriertransform of the oblique sinograms. See Defrise, Michel, A FactorizationMethod for the 3D x-ray Transform, Inverse Problems, Vol. 11, pp. 983-94(1995). The present invention extends the derivation of the Fourierre-binning algorithm to TOF data. We assume that each sinogramcharacterized by z, δ and t is measured over the complete rangesε[−R_(FOV),R_(FOV)] and φε[0, 2π). Scanners usually assemble the datainto sinograms with an angular range φε[0, π). Therefore, a 360°sinogram can be built by assembling two 180° sinograms, using thesymmetryp _(t)(s,φ+π,z,δ)=p _(−t)(−s,φ,z,−δ).  (10)This relation holds when the TOF profile is even, h(t, l)=h(−t,−l), aswill be assumed hereinafter. It will be noted that in contrast with thenon-TOF case, two assembled 360° sinograms of opposite δ are notequivalent, and therefore p_(t)(s, φ, z, δ) and p_(t)(s, φ, z,−δ) yieldtwo independent contributions to the re-binned data. To derive theFourier re-binning approximation, we consider the 2D Fourier transformof one TOF sinogram,

$\begin{matrix}{{P_{t}\left( {\omega,k,z,\delta} \right)} = {\int_{- R_{FOV}}^{R_{FOV}}\;{{\mathbb{d}s}{\int_{0}^{2\pi}\ {{\mathbb{d}{{\phi exp}\left( {{{- {\mathbb{i}}}\;\omega\; s} - {{\mathbb{i}}\; k\;\phi}} \right)}}{p_{t}\left( {s,\phi,z,\delta} \right)}}}}}} & (11)\end{matrix}$where ωε

is the radial frequency conjugate to s, and kεZ is the integer azimuthalfrequency conjugate to φ. Replacing p_(t) by its definition (3) andtransforming integration variables (s, l)→(x, y), we get

$\begin{matrix}{{P_{t}\left( {\omega,k,z,\delta} \right)} = {\int_{{\mathbb{R}}^{2}}^{\;}{{dx}\ {\mathbb{d}y}{\int_{0}^{2\pi}\ {{\mathbb{d}{{\phi exp}\left( {- {{\mathbb{i}\Phi}\left( {x,y,\omega,k,\phi} \right)}} \right)}} \times {f\left( {x,y,{z + {\delta\left( {{{- x}\;\sin\;\phi} + {y\;\cos\;\phi}} \right)}}} \right)}{h\left( {t,{{{- x}\;\sin\;\phi} + {y\;\cos\;\phi}}} \right)}}}}}} & (12)\end{matrix}$where the phase of the complex exponential isφ(x,y,ω,k,φ)=ω(x cos φ+y sin φ)+kφ  (13)

The stationary phase approximation is based on the observation that thephase φ varies rapidly as a function of the integration variable φ whenthe frequencies ω and k are large. The complex exponential exp(−iφ) in(12) is then a rapidly oscillating sinusoid-like function. If f and hare sufficiently smooth functions of φ, this behavior results in acancellation of the φ integral over each period of the oscillation. Themajor contribution to the φ integral then comes from the values of φwhere the phase is extremum (stationary), because the oscillations areminimized in the vicinity of these extrema. The extrema of the phase arethe solutions of

$\begin{matrix}{\frac{\partial{\Phi\left( {x,y,\omega,k,\phi} \right)}}{\partial\phi} = {{{\omega\left( {{{- x}\;\sin\;\phi} + {y\;\cos\;\phi}} \right)} + k} = 0}} & (14)\end{matrix}$or−x sin φ+y cos φ=−k/ω.  (15)Thus, even though the values of φ for which the phase Φ is extremum dodepend on x and y, they always correspond to the same position l alongthe LORs. To a first approximation, we can therefore replace l=−x sinφ+y cos φ by −k/ω in the third argument off in equation (12):

$\begin{matrix}{{P_{t}\left( {\omega,k,z,\delta} \right)} \simeq {\int_{{\mathbb{R}}^{2}}^{\;}{{dx}\ {\mathbb{d}y}{\int_{0}^{2\pi}\ {{\mathbb{d}{{\phi exp}\left( {{- {\mathbb{i}}}\;{\Phi\left( {x,y,\omega,k,\phi} \right)}} \right)}} \times {f\left( {x,y,{z - \frac{k\;\delta}{\omega}}} \right)}{{h\left( {t,{{{- x}\;\sin\;\phi} + {y\;\cos\;\phi}}} \right)}.}}}}}} & (16)\end{matrix}$

Noting that the TOF profile is independent of z and δ, and comparingwith equation (12) at δ=0, we obtain the Fourier re-binning equation forTOF PET data,

$\begin{matrix}{{P_{t}\left( {\omega,k,z,\delta} \right)} \simeq {{P_{t}\left( {\omega,k,{z - \frac{k\;\delta}{\omega}},0} \right)}.}} & (17)\end{matrix}$This approximate relation is identical to the standard Fourierre-binning for non-TOF data, applied separately to each time bin t(after merging with the opposite bin −t as discussed above). Re-binningis then achieved by averaging for each slice z₀ the estimates ofP_(reb,t)(ω, k, z₀)=P_(t)(ω, k, z₀, 0) provided by equation (17) for allavailable values of δ.The stationary phase approximation holds asymptotically for |ω|→∞ and|k|→∞. At low frequencies, the accuracy of equation (17) breaks down.Therefore, as in the standard Fourier re-binning algorithm, the lowfrequencies are re-binned using SSRB,P _(t)(ω,k,z,δ)≈P _(t)(ω,k,z,0)  (18)incorporating only the smallest values of δ, for which thisapproximation is accurate. The efficiency of TOF-FORE stems from theempirical observation that the relation defined by equation (17) isaccurate even for low frequencies, so that equation (18) needs only beused for a few frequency samples around the DC term ω=k=0.

Next, consistency conditions are derived for 3D TOF data with a Gaussianmodel of the TOF profile,h(t,l)=exp(−(t−l)²/2σ²)  (19)where σ, the standard deviation of the TOF profile, is related to thefull-width at half-maximum T_(FWHM) of the time difference measurementby

$\begin{matrix}{\sigma = \frac{{cT}_{FWHM}}{4\sqrt{2\log\; 2}}} & (20)\end{matrix}$Inserting the Gaussian model (19) into equation (3), it can be shownthat any function p_(t)(s, φ, z, δ) that can be represented by equation(3) for some twice continuously differentiable function f(x, y, z), mustbe a solution of the two partial differential equations,

$\begin{matrix}{{{\frac{\partial^{2}p_{t}}{{\partial z}\;{\partial\phi}} + \frac{\partial^{2}p_{t}}{{\partial s}{\partial\delta}}} = {{{- s}\;\delta\frac{\partial^{2}p_{t}}{\partial z^{2}}} - {\frac{st}{\sigma^{2}}\frac{\partial p_{t}}{\partial z}} + {\frac{s}{\sigma^{2}}\frac{\partial p_{t}}{\partial\delta}}}}{and}} & (21) \\{{{{- t}\frac{\partial p_{t}}{\partial z}} + \frac{\partial p_{t}}{\partial\delta}} = {\sigma^{2}{\frac{\partial^{2}p_{t}}{{\partial z}{\partial t}}.}}} & (22)\end{matrix}$When σ→∞, the two last terms on the RHS of equation (21) vanish, andthat equation reduces to John's Equation for non-TOF data (John, F., TheUltrahyperbolic Equation With Four Independent Variables, Duke Math. J.Vol. 4, pp. 300-22 (1938)). For simplicity we will also refer toequation (21) as “John's Equation.” In the opposite limit where σ→0,both equation (21) (if s≠0) and equation (22) reduce to

$\begin{matrix}{{{{- t}\frac{\partial p_{t}}{\partial z}} + \frac{\partial p_{t}}{\partial\delta}} = 0} & (23)\end{matrix}$which is equivalent to the TOF-SSRB relation in equation (8).

The existence of two independent consistency conditions (21) and (22),instead of only one for non-TOF 3D data, can be understood by notingthat four parameters are required to parameterize non-TOF data, whereasthe TOF data depend on five parameters (t, s, φ, z, δ), and f(x, y, z)in both cases only depends on three parameters. Two approaches to exactre-binning will now be derived from the two consistency conditions.

We first focus on John's Equation (21), which can be applied separatelyto each time bin t. Taking the 2D Fourier transform of equation (21)with respect to s and φ, and using the equivalence between thederivative of a function with respect to s (or φ) and the multiplicationof its Fourier transform by iω (or ik), leads to the following equationfor the function Pt(ω, k, z, δ) defined in (11):

$\begin{matrix}{{{k\frac{\partial P_{t}}{\partial z}} + {\omega\frac{\partial P_{t}}{\partial\delta}}} = {{{- \delta}\frac{\partial^{3}P_{t}}{{\partial\omega}{\partial z^{2}}}} - {\frac{t}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial z}}} + {\frac{1}{\sigma^{2}}{\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial\delta}}.}}}} & (24)\end{matrix}$An exact re-binning algorithm similar to the FOREJ algorithm describedin “Fast Rebinning Algorithm for 3D PET Using John's Equation,” Defrise,M. and Liu, X., Inverse Problems, Vol. 15, pp. 1047-65 (1999), isobtained by considering a fixed (ω≠0, k), and by noting that the LHS ofequation (24) is the directional derivative of P_(t) along the vector(k, ω) in the plane (z, δ). In this plane a line segment z=z₀+kδ/ω,0≦δ≦δ₁, is defined, which links a point (z₁=z₀+(k/ω)δ₁, δ₁≠0),corresponding to some measured oblique sinogram, to a point (z₀, 0)corresponding to the 2D sinogram of slice z₀. Along this line segment,equation (24) is written as

$\begin{matrix}{\frac{\mathbb{d}{P_{t}\left( {\omega,k,{z_{0} + {k\;{\delta/\omega}}},\delta} \right)}}{\mathbb{d}\delta} = {{- \frac{1}{\omega}}\left( {{\delta\frac{\partial^{3}P_{t}}{{\partial\omega}{\partial z^{2}}}} + {\frac{t}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial z}}} - {\frac{1}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial\delta}}}} \right){\left( {\omega,k,{z_{0} + \frac{k\;\delta}{\omega}},\delta} \right).}}} & (25)\end{matrix}$Integrating between δ=0 and δ=δ₁ leads to the exact re-binning equation

$\begin{matrix}\begin{matrix}{{P_{t}\left( {\omega,k,z_{0},0} \right)} = {{P_{t}\left( {\omega,k,z_{1},\delta_{1}} \right)} -}} \\{\int_{0}^{\delta_{1}}\mspace{7mu}{{\mathbb{d}\delta}\frac{\mathbb{d}{P_{t}\left( {\omega,k,{z_{0} + {k\;{\delta/\omega}}},\delta} \right)}}{\mathbb{d}\delta}}} \\{= {{P_{t}\left( {\omega,k,{z_{0} + {k\;{\delta_{1}/\omega}}},\delta_{1}} \right)} +}} \\{\frac{1}{\omega}{\int_{0}^{\delta_{1}}\mspace{7mu}{\mathbb{d}{\delta\left( {{\delta\frac{\partial^{3}P_{t}}{{\partial\omega}{\partial z^{2}}}} +} \right.}}}} \\\left. {{\frac{t}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial z}}} - {\frac{1}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial\delta}}}} \right) \\{\left( {\omega,k,{z_{0} + \frac{k\;\delta}{\omega}},\delta} \right).}\end{matrix} & (26)\end{matrix}$As in the previous section, re-binning is achieved by averaging for eachslice z₀ the estimates of P_(reb,t)(ω, k, z₀)=P_(t)(ω, k, z₀, 0)provided by (26) for all available values of δ₁.

If the axial aperture δ is sufficiently small, and the standarddeviation σ of the TOF profile is sufficiently large, the RHS ofequation (25) can be neglected. At this, approximation the Fouriertransformed sinogram P_(t) is constant along the line z=z₀+kδ/ω, theintegral on the RHS of equation (26) disappears and equation reduces toequation (17). This alternative derivation of TOF-FORE shows that theapproximation error has contributions of order O(δ²/ω) and of orderO(δ/σ²ω).

An alternative approach to exact re-binning is based on the secondconsistency condition (22). Consider a fixed (s, φ) and define in theplane (z, δ) a line segment z=z₀−tδ, 0≦δ≦δ₁, which links a point(z₁=z₀−tδ₁, δ₁≠0), corresponding to some measured oblique sinogram, to apoint (z₀, 0) corresponding to the 2D sinogram of slice z₀. Thedirectional derivative of p_(t) along this line is

$\begin{matrix}\begin{matrix}{\frac{\mathbb{d}{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta}},\delta} \right)}}{\mathbb{d}\delta} = {\left( {{{- t}\frac{\partial p_{t}}{\partial z}} + \frac{\partial p_{t}}{\partial\delta}} \right)\left( {s,\phi,{z_{0} - {t\;\delta}},\delta} \right)}} \\{= {\sigma^{2}\frac{\partial^{2}{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta}},\delta} \right)}}{{\partial z}{\partial t}}}}\end{matrix} & (27)\end{matrix}$where we have used the consistency condition (22); Integrating betweenδ=0 and δ₁ leads to

$\begin{matrix}{{p_{t}\left( {s,\phi,z_{0},0} \right)} = {{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta_{1}}},\delta_{1}} \right)} - {\sigma^{2}{\int_{0}^{\delta_{1}}\mspace{7mu}{{\mathbb{d}\delta}{\frac{\partial^{2}{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta}},\delta} \right)}}{{\partial z}{\partial t}}.}}}}}} & (28)\end{matrix}$

A few differences with the first exact re-binning, equation (26), areworth noting:

-   -   Re-binning is applied directly to the sinograms, without        requiring a 2D Fourier transform.    -   The TOF sampling should be sufficiently fine to allow an        accurate calculation of the partial derivative with respect to        t.

The TOF-SSRB equation (8) is obtained by neglecting the integral on theRHS of equation (28), which is an approximation of order O(σ²δ). Thisshows that the accuracy of TOF-FORE and of TOF-SSRB have oppositebehaviors when the TOF resolution improves. For the relatively poortiming resolution achievable with current detectors, TOF-FORE appearspreferable, as demonstrated by the numerical results below.

TABLE 1 Simulation parameters. Ring radius (mm) 278 Number of rings 32Maximum ring difference used 29 Span 1 Number of Oblique sinograms 1018tan θ_(max) (degree) 14 Ring spacing (mm) 4.8 Radial sampling (mm) 1.2Number of radial s samples 256 Number of angular φ samples 256 Digitizedphantom 256 × 256 × 63 voxel size (mm) 1.2 × 1.2 × 2.4The TOF-FORE algorithm has been evaluated using simulated data for amulti-ring scannerdescribed by the parameters in table 1. Experimentally, TOF data (with atime resolutionof around 1.2 ns) are sorted in multiple time bins (e.g. 9 bins of 500ps) to cover the wholeFOV (4.5 ns=675 mm). In this simulation, we examine TOF data with aGaussian profile ofFWHM 75 mm or 120 mm, centered at t=50 mm. This corresponds to a timeresolution of500 ps or 800 ps at Δτ=300 ps. Noise free 3D TOF sinograms weresimulated according toequation (1), by forward projecting the digitized image of two phantoms(without attenuationor scatter). The modified 3D TOF data (3) were then obtained using theapproximation (4).The oblique sinograms were re-binned with SSRB, TOF-SSRB and TOF-FORE.For TOF-FOREthe low-frequency samples (i_(ω), k) satisfying |i_(ω)|≦i_(max)=6 and|k|≦k_(max)=6 were re-binned using SSRB with a maximum ring differencerd_(max)=4. The re-binned 2D sinograms were compared to the ‘exact’ 2Dsinograms p_(exact,t)(s, φ, z) obtained by forward projecting thedigitized phantom with the same TOF profile. For each slice, anormalized root-mean-square (RMSE) difference between the re-binnedsinogram and the ‘exact’ sinogram was calculated as

$\begin{matrix}{{{RMSE}\lbrack z\rbrack} = {\left( \frac{\int\;{{\mathbb{d}s}{\int\;{\mathbb{d}{\phi\left( {{p_{{reb},t}\left( {s,\phi,z} \right)} - {p_{{exact},t}\left( {s,\phi,z} \right)}} \right)}^{2}}}}}{\max_{z}{\int\;{{\mathbb{d}s}{\int\;{{\mathbb{d}\phi}\;{p_{{exact},t}\left( {s,\phi,z} \right)}^{2}}}}}} \right)^{1/2}.}} & (29)\end{matrix}$For comparison, non-TOF data of the two phantoms were also simulated andre-binned with the standard FORE algorithm, using the same parametersi_(max), k_(max) and rd_(max) for the low frequency region. The RMSE wasthen calculated according to (29) (normalization is important becausethe sinogram values are smaller for the TOF data, due to multiplicationby the TOF profile (19)).

The first phantom is an ellipsoid with uniformly distributed activity,centered at (x, y, z)=(56,−76, 75) mm with semi-axis (20, 30, 3) mm. Thecenter of this object is not aligned axially with the center of areconstructed slice. FIG. 4 shows the RMSE difference between there-binned and the exact sinograms, as a function of the axial slice. Thesecond phantom is a combination of 20 cylinders of various sizes andintensities, inside a larger cylinder of radius 120 mm (see FIG. 5). Forthis phantom, the RMSE was normalized with the root-mean-square of eachslice (i.e. using equation (29) without the max_(z)). FIG. 6 shows theRMSE for the various algorithms with FWHM=120 mm, plotted versus thenormalized RMSE for the non-TOF data re-binned with FORE. On theaverage, the error with TOF-FORE is a factor 1.38 larger (slope of theregression line) than the non-TOF reference (solid line in FIG. 6). Theerror does not exceed 5% and is significantly smaller than with the twoSSRB algorithms. As with the ellipsoid phantom, using the TOFinformation in SSRB significantly improves the accuracy, as shown inTable 2 below.

TABLE 2 Equations of the regression lines in FIGS. 6 and 7 FWHM = 75 mmFWHM = 120 mm SSRB y = 3.23x + 0.012 y = 2.62x + 0.009 TOF-SSRB y =2.46x + 0.002 y = 2.26x + 0.006 TOF-FORE y = 1.71x + 0.000 y = 1.38x +0.000

As expected from the analysis of John's equation above, the accuracy ofTOF-FORE decreases when the TOF resolution improves. For instance, for aFWHM of 75 mm (FIG. 7), the error with TOF-FORE is a factor 1.71 largerthan the non-TOF reference, instead of 1.38 with a FWHM of 120 mm. FIGS.8( a)-(d) and 9(a)-(d) illustrate the structure of the re-binning errorfor this phantom.

In FIGS. 8( a)-(d), the re-binned sinogram i_(z)=32 of the cylindricalphantom is obtained from TOF data with FWHM=120 mm. FIG. 8( a) shows the“exact” sinogram. The absolute error for the TOF-SSRB re-binning, scaledto 10% of the maximum of the “exact” sinogram, is shown in FIG. 8( b).The absolute error for TOF-FORE, scaled to 10%, is shown in FIG. 8( c).The same absolute error for TOF-FORE, scaled to 2% to elicit low-levelartifacts, is shown in FIG. 8( d).

FIGS. 9( a)-(d) illustrate an axial section of the stack of re-binnedsinograms for the cylindrical phantom with TOF data (FWHM=120 mm),respectively for the “exact” sinogram, the TOF-FORE re-binding, TOF-SSRBre-binning, and SSRB re-binning. The section shown corresponds to s=−84mm. The vertical axis is the axial slice coordinate z and the horizontalaxis is the angle φ. Grey scale range [0.40M, M], where M is the maximumvalue in the ‘exact’ axial section.

In summation, measuring and exploiting time-of-flight information in 3DPET has considerable potential, especially for whole-body imaging oflarge patients, an application where improvement would have an importantclinical impact. The present invention provides a way to accelerate thereconstruction of 3D TOF data through Fourier re-binning. The TOF-FOREalgorithm is provided to generate a good approximation of “exact” 2D TOFsinograms and shown to be superior to TOF-SSRB when the timingresolution is larger than 500 ps and the axial aperture does not exceed15°. Simulation results demonstrate that TOF-FORE can be a valuablemethod for PET scanners with TOF capability: the timing resolutionachievable with current detector technology is larger than 1 ns, and inthis case the accuracy of TOF-FORE is comparable to the accuracy of theFORE algorithm.

The present invention further provides two partial differentialequations that must be satisfied by consistent TOF-PET data. These twoequations lead to two different exact re-binning algorithms, whichprovide insight into the accuracy of the TOF-FORE and TOF-SSRBalgorithms, both of which were shown to be first-order approximations tothe exact re-binning.

The five dimensional data space in 3D TOF-PET results in a very richstructure which was only partially elicited and exploited here by thederivation of the two consistency conditions. It will be apparent tothose skilled in the art from the present disclosure that more efficientre-binning algorithms, which would better combine the measured TOFinformation with the ‘virtual’ TOF information provided by thestationary phase approximation, may be achieved by further analysis ofthe 3D TOF-PET data structure.

It should be appreciated by those having ordinary skill in the art thatwhile the present invention has been illustrated and described in whatis deemed to be the preferred embodiments, various changes andmodifications may be made to the invention without departing from thespirit and scope of the invention. Therefore, it should be understoodthat the present invention is not limited to the particular embodimentsdisclosed herein.

1. A method for reconstructing a nuclear medical image from TOF-PETimaging data, comprising the steps of: obtaining three dimensionalTOF-PET data from a PET scanner; storing said TOF-PET data in aplurality of TOF time bins, each bin corresponding to a time differencevalue of photons arriving at opposite detectors of said PET scanner;estimating, for each time bin, a two dimensional TOF sinogram of eachtransaxial slice z within an axial field of view of said scanner, byaveraging independent estimates of redundant three dimensional TOF-PETdata; and reconstructing an image from said estimated two dimensionalTOF sinograms using a two dimensional TOF reconstruction algorithm. 2.The method of claim 1, wherein said three dimensional TOF-PET data isrepresented by the following equation:${p_{t}^{m}\left( {s,\phi,z,\delta} \right)} = {\sqrt{1 + \delta^{2}}{\int_{- \infty}^{\infty}\mspace{7mu}{{\mathbb{d}{{lf}\left( {{{s\;\cos\;\phi} - {l\;\sin\;\phi}},{{s\;\sin\;\phi} + {l\;\cos\;\phi}},{z + {l\;\delta}}} \right)}}{{h\left( {t,{l\sqrt{1 + \delta^{2}}}} \right)}.}}}}$3. The method of claim 2, wherein the step of estimating comprises thestep of modifying said three dimensional TOF-PET data top_(t)(s, ϕ, z, δ) = ∫_(−∞)^(∞) 𝕕lf(s cos  ϕ − l sin  ϕ, s sin  ϕ + l cos  ϕ, z + l δ)h(t, l).4. The method of claim 3, wherein said modified data is estimated asp _(reb,t)(s,φ,z+tδ)=p _(t)(s,φ,z,δ).
 5. The method of claim 4, whereinsaid two dimensional TOF sinogram is obtained by solving for theequation${p_{{reb},t}\left( {s,\phi,z} \right)} \simeq {\frac{1}{2{\delta_{\max}(z)}}{\int_{- {\delta_{\max}{(z)}}}^{\delta_{\max}{(z)}}{{\mathbb{d}\delta}\;{{p_{t}\left( {s,\phi,z,\delta} \right)}.}}}}$6. The method of claim 5, wherein said two dimensional TOF sinogram isobtained by solving for the equation${P_{t}\left( {\omega,k,z,\delta} \right)} \simeq {{P_{t}\left( {\omega,k,{z - \frac{k\;\delta}{\omega}},0} \right)}.}$7. The method of claim 5, wherein said two dimensional TOF sinogram isobtained by solving for the equation $\begin{matrix}{{P_{t}\left( {\omega,k,z_{0},0} \right)} = {{P_{t}\left( {\omega,k,z_{1},\delta_{1}} \right)} - {\int_{0}^{\delta_{1}}{{\mathbb{d}\delta}\frac{\mathbb{d}{P_{t}\left( {\omega,k,{z_{0} + {k\;{\delta/\omega}}},\delta} \right)}}{\mathbb{d}\delta}}}}} \\{= {{P_{t}\left( {\omega,k,{z_{0} + {k\;{\delta_{1}/\omega}}},\delta_{1}} \right)} +}} \\{\frac{1}{\omega}{\int_{0}^{\delta_{1}}{\mathbb{d}{\delta\left( {{\delta\frac{\partial^{3}P_{t}}{{\partial\omega}{\partial z^{2}}}} + {\frac{t}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial z}}} - {\frac{1}{\sigma^{2}}\frac{\partial^{2}P_{t}}{{\partial\omega}{\partial\delta}}}} \right)}}}} \\{\left( {\omega,k,{z_{0} + \frac{k\;\delta}{\omega}},\delta} \right).}\end{matrix}$
 8. The method of claim 5, wherein said two dimensional TOFsinogram is obtained by solving for the equation${p_{t}\left( {s,\phi,z_{0},0} \right)} = {{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta_{1}}},\delta_{1}} \right)} - {\sigma^{2}{\int_{0}^{\delta_{1}}{{\mathbb{d}\delta}{\frac{\partial^{2}{p_{t}\left( {s,\phi,{z_{0} - {t\;\delta}},\delta} \right)}}{{\partial z}{\partial t}}.}}}}}$